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In this article we study the classical nearest-neighbour spin-ice model (nnSI) by means of Monte 
Carlo simulations, using the Wang-Landau algorithm. The nnSI describes several of the salient 
features of the spin-ice materials. Despite its simplicity it exhibits a remarkably rich behaviour. 

The model has been studied using a variety of techniques, thus it serves as an ideal benchmark to 
test the capabilities of the Wang Landau algorithm in magnetically frustrated systems. We study in 
detail the residual entropy of the nnSI and, by introducing an applied magnetic held in two different 
crystallographic directions ([111] and [100],) we explore the physics of the kagome-ice phase, the 
transition to full polarisation, and the three dimensional Kasteleyn transition. In the latter case, 
we discuss how additional constraints can be added to the Hamiltonian, by taking into account 
a selective choice of states in the partition function and, then, show how this choice leads to the 
realization of the ideal Kasteleyn transition in the system. 


I. INTRODUCTION 

In magnetism, when a spin cannot fully minimise its in¬ 
teractions with its neighbours, the system is called “frus¬ 
trated” . This situation can arise under a number of 
different circumstances, such as bond disorder, further 
neighbour interactions and lattice geometry. Geometri¬ 
cally frustrated magnets are the cleanest and better con¬ 
trolled in experimental systems. Frustration precludes 
the formation of simple ordered ground states, rather, it 
typically leads to a degenerate manifold of ground states, 
scaling with the system size, therefore, exhibiting ex¬ 
tensive zero-point entropy. Unsurprisingly, these ground 
states are unstable to any small perturbation. Frustrated 
systems exhibit a rich variety of behaviour, including or¬ 
der by disorder, fractionalisation and magnetic analogues 
of solids, liquids, glasses, ice, quantum liquids and bose 
condensation. They represent ideal model systems for the 
study of generic concepts relevant to collective phenom¬ 
ena, where simple classical Hamiltonians can give rise to a 
wealth of different phenomenai^— . Given that analytical 
treatment cannot be performed, frustrated magnets are 
very well suited for numerical modelling, and constitute 
ideal test-grounds for powerful numerical techniques. 

Among the recent developments in classical simula¬ 
tion techniques, the Wang-Landau algorithm (WLA) has 
stood out as a powerful tool for the determination of 
phase diagrams^. The aim of this technique is to give 
a very accurate estimate of the density of states of the 
system (DoS), from which the thermodynamic behaviour 
can be calculated using canonical ensemble statistical me¬ 
chanics. There are two particular points of this technique 
that makes it attractive for the study of frustrated sys¬ 


tems: first, by simply providing an estimate of the DoS 
it gives direct information about the degeneracy of the 
ground state and its low temperature excitations; second, 
in its algorithmic construction it has a built in mechanism 
to avoid the trapping of the system into local minima 
in the energy landscape. The existence of many local 
minima in frustrated magnets is one of the main rea¬ 
sons of the under-performance of the usual Monte Carlo 
techniques. Recently, the WLA has been successfully ap¬ 
plied to specific problems in frustrated systems and dimer 
models, such as the formation of magnetisation plateaus 
in the Ising Shastry-Sutherland model^ or to determine 
the order of a transition in the Heisenberg stacked trian¬ 
gular antiferromagnelji and in dimer models with next- 
nearest-neighbor interactions^. In the present work we 
explore thoroughly a simple classical frustrated model, 
the nearest-neighbour spin-ice model (nnSP), by means 
of the WLA. This model, despite its apparent simplicity 
—Ising spins on a pyrochlore lattice interacting ferromag- 
netically at nearest neighbours— exhibits an extremely 
rich behaviour. The degeracy of the zero field ground 
state grows exponentially with the system size and with 
aufbau rules analogous to those of protons in water ice. 
By applying an external magnetic field it is possible to 
find a metamagnetic transition or to tune into an effec¬ 
tive two dimensional model (the kagome-ice, also with 
extensive residual entropy) and to find two- and three- 
dimensional Kasteleyn transitions. The physics of this 
model system has been well studied using a variety of 
techniques^rJ^, and thus serves as an ideal benchmark 
to test the capabilities of the WL algorithm in this type 
of systems. Beyond that, the WLA provides new ther¬ 
modynamic results that have not been obtained by other 
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methods, such as the free energy as a function of the order 
parameter, which, to our knowledge, has been calculated 
for a Kasteleyn transition for a first time. 

The structure of this article is as follows. In the next 
section we give a brief discussion of the model and the 
simulation technique for completeness. In section III, 
we study the model at zero applied field, calculate the 
residual entropy and compare it with different existing 
estimates using the WLA technique. In the same section, 
we study the system under field applied along [III], and 
explore the entropy of the kagome-ice state, as well as the 
entropy peak that rises when this state is destroyed by 
increasing the temperature. Finally, we study the case 
of field applied along [100], and do a characterisation of 
the three dimensional Kasteleyn transition into the fully 
polarised state. 

II. METHODS 
A. Wang Landau Algorithm 

Recently, F. Wang and D. P. Landau introduced an 
algorithm to estimate the density of states of a sys¬ 
tem by performing a random walk in energy space^i^. 
This method is closely related to ‘umbrella sampling’ 
techniques!^ and multicanonical Monte Carloi^. The al¬ 
gorithm provides a very good estimate of the DoS of the 
system over a bounded region of the energy spectrum. 
The DoS can be calculated as a function of the energy, 
if one works in the canonical ensemble, but also as a 
function of other variables like pressure or magnetisa¬ 
tion if one is interested in other ensembles such as the 
isothermal-isobaric or its magnetic equivalent (as we will 
use in sections IlIIBI and llll Cp . In this section we describe 
the procedure for one variable, which is then straightfor¬ 
wardly extended to the case of several variables. 

The algorithm requires the knowledge of the Hamilto¬ 
nian of the system and a method to sample configurations 
- in our case it is simply random single spin flips. The 
starting point is an arbitrary initial configuration with 
energy Eq. Initially, the DoS is taken as homogeneous: 
g{E) = 1 for all E. One step of the calculation consists 
then in choosing new random configuration, calculating 
its energy Ei , and accepting it or discarding it with prob¬ 
ability 

p[Eq^ Ei) = TEm{l,j^^). ( 1 ) 

At each step, the histogram H{E) and the DoS g{E) of 
the final configuration are modified according to g{E) —>■ 
giE)f and H{E) —>■ H{E) -|- 1. Initially the modification 
factor / for g{E) is taken from a larger value (usually, 
/o = e^) which is then reduced as the algorithm pro¬ 
gresses. The rule by which the modification factor fi is 
reduced is an important choice that conditions the ac¬ 
curacy of the DoS and speed at which the process con¬ 
verges. In principle, one can choose any function that 


tends monotonically to one, and stop the process once / 
reaches a given value. In the original work, Wang and 
Landau propose the rule fk+i = '/Jk, but other choices 
are possible which give faster convergence. In particular, 
Belardinelli and Pereyra (BP)!^ proposed and alternative 
method with improved convergence times, in which the 
modification factor is eventually scaled as the inverse of 
the Monte Carlo time, t. For this work we have adopted 
the BP algorithm. 

In our case, within one Monte Carlo step this proce¬ 
dure was repeated with single spin flips until the spin of 
each site had been chosen to change at least once on av¬ 
erage. In practice, the quantity lng(i?) —Ing(if) -I- In/ 
is used for simplicity (hence the choice /o = e). Notice 
that this procedure guarantees that the random explo¬ 
ration in phase space will not jam at local minima: each 
time the transition to the new state is not accepted, the 
DoS of the initial state is increased and thus the prob¬ 
ability of accepting any future transition, which will be 
proportional to g{Eo), is also increased. The detailed 
balance condition is satisfied to within In / accuracy. 

During this random walk through phase space, the his¬ 
togram is accumulated and checked periodically. When 
H{E) becomes flat, the histogram is reset, and the next 
random walk begins with a finer modification factor /i. 
The criterion of flatness varies according the size and 
complexity of the system. Usually, the criterion used for 
flatness is that every entry in H{E) is not smaller than 
a percentage of the average histogram for all E. 

The final result is a relative DoS. To calculate the ab¬ 
solute values one needs some additional information, e.g. 
a known point in the DoS - in our case the high temper¬ 
ature value of the DoS tends to 2, or the knowledge of 
the integral of the DoS -in our case 2^. Once the DoS 
of the system is known, it is straightforward to calcu¬ 
late the thermodynamic quantities in the canonical en¬ 
semble: for example, from Z — ^the 
free energy can be written as F{T) = —kT log Z, the 
internal energy U{T) = ^ ^he entropy 

S{T) = and the specific heat, using the usual 

linear fluctuation relation, C{T) = _ 

B. Spin ice 

The nnSI Hamiltonian under an external magnetic field 
reads: 

H = Jeff ^ S,.S,(2) 
<ij> i 

where the Si are Ising spins situated at the corners of a 
pyrocholore lattice (see inset of Fig. [1]), H is the exter¬ 
nal magnetic field, g the gyromagnetic ratio and JeS is 
the effective exchange interaction which is taken to be 
positive. 

The model is applicable to a certain class of materials, 
of which the most notable examples are Dy 2 Ti 207 and 
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Ho 2 Ti 207 - In these materials, the magnetic ions sit at 
the corners of a pyrochlore lattice and are constrained by 
the crystalline field to point along the local (111) quanti¬ 
zation axes {i. e. to or from the centre of the tetrahedra). 

The nnSI model provides a very good description of these 
materials between 0.2K and lOK— . The main difference 
arises from the fact that in the materials, in addition to tq' 
the exchange interaction - which is antiferromagnetic- ^ 
there is a large long range dipolar interaction. If the lat¬ 
ter is truncated beyond the nearest-neighbour spins, Jgff 
is effectively ferromagnetioi^. 

The ground state of the nnSI scales exponentially with 
system size and obeys the local construction rule that 
within any tetrahedra, two of its spins should point in¬ 
wards, and two outward o^b^^ . This rule is called the 
“ice-rule” given its analogy to the Bernal and Fowler rules 
for protons in water-ice^^ and is the origin of the epithet 
“spin-ice” given to these models. The exponential degen¬ 
eracy of this ground state leads to an extensive residual 
entropy at zero temperature. 

We analysed the nnSI model using the WLA with both 
the original implementation and the one proposed by Be- 
lardinelli and Pereyra. In order to reach larger lattices 
we also performed simulations dividing the energy range 
in multiple regions. To normalize the DoS, we generally 
used the condition that there are only two states with 
maximum energy (the “all-in” and “all-out” configura¬ 
tions). The comparison between several runs performed 
with different seeds, or different normalisation, allowed us 
to estimate the errors on the residual entropy for differ¬ 
ent sizes. We explored the configuration space through 
random single spin-flip moves. In addition, we used a 
conventional cubic cell for the pyrochlore lattice, which 
contains 16 spins, and simulated systems with Lx L cells, 
with L ranging from 1 to 8 (16 to 8192 spins). The DoS 
was estimated as a function of energy with the modifica¬ 
tion factor changing from /q = to /final = exp(10“®). 

For the results of the section ITlI Al we used the single¬ 
variable WLA, while for sections IIII Bl and IIII Cl we ac¬ 
cumulated the DoS as a function of both the energy and 
the magnetisation in the direction of the applied mag¬ 
netic field. In this latter case, the thermodynamic quan¬ 
tities are calculated using a partition function which is 
summed in the energy and in a variable conjugate to the 
magnetic field, in this way the derivatives of Z provide 
information about the average magnetisation. We have 
chosen J = 1.11 AT in order to match the effective nn 
exchange constant for Dy 2 Ti 207 . 


III. RESULTS 

A. Spin-ice with no applied external field: residual 
entropy 

As a first step we analysed the nnSI model with no 
applied external field {H = 0), in particular, we were 
interested in the behaviour of the specific heat and the 



E/k (K/Dy ion) 

FIG. 1. Natural logarithm of the density of states of the nnSI 
model for a system size L = 4. Its shape is characteristically 
asymmetric with high degeneracy of the ground state. The 
inset shows an schematic view of the pyrochlore lattice. 

entropy. In Fig.[T]it is shown the logarithm of the DoS of 
the nnSI model for a system size L = 4. Its shape, in con¬ 
trast with the highly symmetrical DoS of the Ising model 
(see e. g. d), is characteristically asymmetric, starting 
from a high value at the lowest energy, a direct measure 
of the high degeneracy of the ground state. 

In the case of the specific heat, the expectation is that 
as the system is cooled down, there will be an onset of 
correlations as the temperature T approaches J/k, that 
will be evidenced by a Schottky-like peaki^. The spe¬ 
cific heat calculated using the WLA shows a peak close 
to T = IK (green crosses in figure [2]), in accordance 
with those calculated by other techniques^ and with the 
high temperature experimental results in spin-ice mate¬ 
rials (T > 0.6K)^. This same figure shows the tempera¬ 
ture dependence of the entropy per mole of the system. 
The high temperature value (not shown in the figure) is 
i?ln2 Ri 5.76 J/molK, indicating that the system is be¬ 
having as an uncorrelated paramagnet. As the tempera¬ 
ture is lowered the entropy is reduced until it reaches a 
residual value S'o close to 1.7 J/molK. 

The residual entropy is a characteristic feature of ice 
models; in real spin-ice materials, such as Dy 2 Ti 207 or 
Ho 2 Ti 207 it is expected that the degeneracy of the spin- 
ice manifold will be lifted by additional interactions — 
chiefly the dipolar interaction— and that the system at 
T = 0 will be ordere d^^i^^ . In the nnSI model, however, 
the aufbau rules are strictly those of Bernal and Fowler, 
and one expects to find an exponentially degenerate state 
with the same extensive residual entropy of the three 
dimensional ice models. The determination of the value 
of this entropy has a long history, starting with Linus 
Pauling’s famous estimate in 1935^. Sq can be written 
as 


So = kin Wn 


( 3 ) 
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FIG. 2. Specific heat (green) and entropy (red) vs temper¬ 
ature for system size L = 4, calculated from the density of 
states (DoS) with parameters set for Dy 2 Ti 2 O7. The figure 
shows the expected Schottky-like peak in the specific heat, 
marking the onset of spin-ice correlations. The entropy does 
not vanish as the temperature is lowered towards T — 0 but 
shows the residual value characteristic of all ice models. 


where k is Boltzmann’s constant and Wn = is the 
number of microstates that form the ground state of the 
system, with Nt the number of tetrahedra. Pauling’s es¬ 
timate gives W = which translated to the entropy 
per mole gives Sq = i?/21n3/2 Ri 1.68 J/molK, in rea¬ 
sonable agreement with the results of Fig[5J 

Pauling’s estimate neglects correlation effects (in the 
form of loops) and it can be shown that it is a lower 
bound on the true While an exact solution exist 

for two dimensional ice models^, this is not true for three 
dimensions. Currently the best estimate of the entropy is 
that due to J. F. Nagle^^, who, building on work by Still- 
inger and Di Marzio^l, used a series expansion method 
to derive the estimate 


bFNagie = 1.50685(15). (4) 

The most accurate calculation of Sq for the nnSI model 
in the literature comes from the integration of energy 
and magnetisation data obtained by loop Monte Carlo 
simulationsiS and gives Wi = 1.5071(3), very close to 
Nagle’s result. The WLA provides a direct determination 
of the entropy without the need of integrating the spe¬ 
cific heat and specifying additional constants, and is thus 
ideally suited for an accurate determination of the resid¬ 
ual entropy of the nnSI model. A variant of the WLA 
has been successfully applied to determine the residual 
entropy of two simple nearest neighbours ice models in 
three dimensional hexagonal lattices: the six-state H 2 O 
molecule model and the two-state H-bond model^^. A 
similar method is applied to nnSI. We have calculated 
the DoS for lattice sizes L = 1 to 8 which correspond to 
Nt = 8 to 4096 tetrahedra, and from those determined 
W for the ground state. Figure [3] shows W as a function 


of the inverse of the number of tetrahedra I/Nt- We 
have used two different criteria to normalise the DoS: in 
one case we used a known density for a given state (the 
highest energy state) and in the other we used the sum 
of all states (2^ in our system). The differences in W 
given by the different criteria of normalisation, or that 
obtained in independent runs using a different set of ran¬ 
dom numbers, are smaller than the size of the symbols 
in the figure. 



1/N^ 


FIG. 3. The number of microstates of the ground state,IT, as 
a function of the inverse of the number of tetrahedra 1 /At. 
As expected, the residual entropy decreases as the system 
size is increased. Shown for reference are the thermodynamic 
value of W from Pauling’s estimate (blue square) and from 
Nagle’s calculation (red circle). The line is a fit according to 
eq. m The inset shows the full range of the fit (from L = 1 to 
L = 8). 

Fig. [3] shows that, as expected, the residual entropy 
decreases as the system size is increased. In order to 
obtain the thermodynamic value of W we fit this data to 
the form 

W{x) = Woo + ai(^-^'^ . (5) 

From this fit we obtain Woe = 1.50682(9), with oi = 
1.557(9) and 9 = 0.883(3). The sub-linear value of 6 
is an indication of bond correlations in the ground-state 
manifold. Our value for W in the thermodynamic limit 
is perfectly consistent with the results by Nagle (see eq. 
a and by previous calculation s . 


B. Spin-ice: Field along [111] 

A remarkable feature of the Sinn model is that an ex¬ 
ternal magnetic field can tune the system into regions of 
different physics. Two notable cases happen when the 
field is oriented (i) along the crystallographic [111] di¬ 
rection, which will be described in this subsection, and 
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(ii) along [100], which will be described in the following 
subsection. To describe the effects of an applied external 
field in the WLA,on one hand the Zeeman term must be 
included in the Hamiltonian (the second term in eq. [2]) 
and, on the other hand, it is more convenient to calculate 
the DoS as a function of two indices: energy and mag¬ 
netisation M, the latter being the quantity conjugate to 
the field. In this way is it possible to work in the mag¬ 
netic equivalent of the isothermal-isobaric ensemble and 
obtain directly from the doubly summed partition func¬ 
tion the Gibbs free energy G{H, T) and the mean value 
of M. The disadvantage of this approach is that it is 
too expensive in calculation resources and, therefore, the 
calculations are usually constrained on smaller system 
sizes. 



FIG. 4. Logarithm of the DoS as a function of E and M along 
[111] calculated for a system of size L = 3. The figure also 
shows the projection of the DoS over the M — E plane. 

Fig. m shows the logarithm of the DoS as a function 
of E and M along [111] calculated for a system of size 
L = 3. As expected, the DoS is still asymmetric along 
the E direction (cf. Fig. [T]) and symmetric in the M axis. 
The figure also shows the projection of the DoS over the 
M — E plane, which takes the shape of a pentagon. In 
this projection one can see that while the highest energy 
state corresponds to M = 0 the ground state manifold 
comprises a range of magnetisation values; for H || [111] 
this range goes from -3.33 to 3.33 hb which is the highest 
value of the magnetisation that can be obtained along 
[111] without breaking the ice rules. 

In Refs 0 and [nl the case of the external held along 
[111] has been studied for the nnSI model both analyti¬ 
cally and through conventional Monte Carlo simulations. 
Along this crystallographic direction the system can be 
thought of as a stack of alternating kagome and trian¬ 
gular planes. In the triangular planes the projection of 
the spins in the [111] direction is 1, while in the kagome 
planes it is 1/3. The evolution of the system under held 
can be obtained by analysing the behaviour of the mag¬ 
netisation at low temperatures. In Fig. [5] we show M as a 
function of the applied magnetic held as calculated using 
the WLA with two variables (energy and magnetisation) 


for L = 3 and for the parameters of Dy 2 Ti 207 , which is 
in perfect coincidence with all the features found in ref. 
0. At very low helds the magnetisation rises linearly 
with a slope which is proportional to 1/T. In this regime 
the held is merely selecting a subset with non-zero mag¬ 
netisation along [111] from the ground-state manifold. 
The slope is given by the competition between the gain 
in Zeeman energy and the entropy, because the num¬ 
ber of states at a given magnetisation decreases sharply 
as M increases. This is shown in the inset of Fig. [5] 
where the logarithm of the DoS is plotted as a function 
of the magnetisation at a hxed energy {k O.IK above 
the ground state). This shows how the logarithm of the 
number of available states is halved as the magnetisation 
is raised towards the 3.33/rB/Dy ion. When the field is 
increased to a value high enough to overcome the en- 
tropic effects, but low enough that the ice-rules are still 
obeyed, the spins in the triangular planes (with projec¬ 
tion 1) all align with the field, the kagome planes decou¬ 
ple, and the system becomes effectively two-dimensional 
termed as “kagome-ice” (KI). The spin-ice rules in the 
KI still allow for an exponentially degenerate number of 
conhgurations, and it still possesses an extensive residual 
entropy (as shown below). The magnetisation reaches a 
plateau in this held range at m = 3.33/iB per Dyspro¬ 
sium ion, which is easily calculated by taking into ac¬ 
count the ice-rules and the different projections of the 
spins in the kagome and triangular planes. If the held is 
increased further, it eventually overcomes the exchange 
interaction and the system goes through a metamagnetic 
transition (at around 1 Tesla in Fig. O to a fully po¬ 
larised state where all spins maximise their projection 
with the magnetic held (m = 5fj,B/^y ion). In the nnSI 
model this transition is merely a crossover and its width 
is strongly temperature dependent. If the dipolar inter¬ 
actions are added to the model, the transition becomes 
a hrst order—. 

As mentioned before, the ice-rules in the KI phase are 
still under-constraining the system and allow for an ex¬ 
ponentially large number of possible conhgurations. In 
this case it is possible to obtain an exact solution for 
the residual entropy due to a mapping from the KI into 
dimers in a honeycomb lattic o^^i^° . This mapping allows 
the study of the effects of slightly misaligned helds, which 
lead to a two-dimensional Kasteleyn transition (see El), 
as well as the transition to the fully polarised state, which 
can be interpreted as a dimer-monomer transition and is 
expected to be accompanied with a peak in the entropy 
as a function of heldiS. 

In conventional MC simulations, the calculation of the 
held dependent entropy usually includes the integration 
of the specihc heat as a function of temperature, with 
the integration constant for each held point determined 
by the value at an appropriate hxed point. In the case 
of the WLA this calculation is straightforward from the 
known g{E,M) and requires no further input. 

We have calculated the entropy as a function of held 
along [111] for diherent temperatures (see Fig. [Q from 
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Magnetic field (T) 


FIG. 5. Magnetization vs field along the [111] axis for L = 
3 and different values of temperature. The magnetization 
curves present two well defined plateaux. The first one, at 
3.33/rs /Dy ion, correspond to a state where all apical spins 
are aligned with the field (the kagome-ice state). The second 
plateau, at 5^b /Dy ion, corresponds to the state of maximal 
spin projection along [111]. The inset shows the logarithm of 
the density of states as a function of the magnetisation at a 
fixed energy {k O.IK above the ground state). 
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FIG. 6. Entropy vs magnetic field for L = 3 and for differ¬ 
ent values of the temperature. The T = 0 curve shows how 
the entropy is reduced to zero through two steps: a first one 
into the KI state, and a second one (at about 1 tesla) when 
the system becomes fully polarised. At higher temperatures 
the most noticeable feature is a giant entropy peak at the 
polarisation transition (see text). 


the DoS for L = 3. The green curve shows the T = 0 
result: here the entropy jumps from its zero field value 
(the 3D ice-entropy) to the KI value (close to 0.7 J/mol 
K) and at higher fields (around I tesla) jumps down to 
zero as the system becomes fully polarised. The red line 
is the exact value for the KI entropy {Sq = 0.6715 J/mol 
K) as calculated in Refs. [Tl| and [SO- The slight discrep¬ 
ancy between the calculated Sq and the exact result is 
due to finite size effects. It is worth pointing out that 
despite its name and the fact that the coordination num¬ 
ber of the lattice is 4, the KI is not sensu stricto an 
ice-type model and thus its residual entropy is smaller: 
Wki = 1.175, compared to Wan = (4/3)3/2 = 1.539... . 
As the temperature increases, the most notable feature is 
the appearance of a giant peak in the entropy, which at 
low temperature (below 0.7K) is larger than the residual 
entropy. It might seem counterintuitive that the applica¬ 
tion of an ordering field might result in an increase of the 
entropy. This can be explained simply, with the aid of the 
dimer mapping, as coming from the crossing of an exten¬ 
sive number of energy levels, corresponding to different 
numbers of monomer defects, which have macroscopic 
entropies (see M). In real materials though, this feature 
of the nnSI model, which has potential applications for 
magnetocaloric manipulations, is almost completely sup¬ 
pressed by additional interactions. 


C. Spin-ice: Field along [100] 

A case of particular interest arises when an external 
field is applied in the [100] direction. In this case, con¬ 


trary to the previous one, the fully saturated state be¬ 
longs to the zero field ground-state manifold, that is, it 
satisfies the ice-rules. At low temperatures {kT <C Jeff) 
and for any value of the magnetic field, there are no exci¬ 
tations in the form of local violations to the ice-rules. In 
Ref. m and [13 it was shown that in this regime the com¬ 
petition between entropy gain and loss of Zeeman energy 
gives rise to a three-dimensional Kasteleyn transition!^ 
where strings of negative magnetisation proliferate and 
span the whole length of the sample. This transition 
takes place at a field dependent critical temperature given 
byi^ 


y/3k B In 2 


( 6 ) 


below which no string is present. This critical tempera¬ 
ture comes from equating in the free energy the loss in 
Zeeman energy per segment of a string of negative mag¬ 
netisation (2/v^h, due to the spin projection) with the 
term arising from the entropic gain per segment (Tin2). 
Since the line spans the whole sample, an equal number 
of in-pointing and out-pointing spins are flipped and the 
ice-rule is preserved in the whole sample, that is to say, 
the Kasteleyn transition occurs between different spin- 
ice states. In the case of our simulation (where we have 
imposed periodic boundary conditions) these lines take 
the shape of non-contractible closed loops in the torus. 
Notice that the characteristic energy of this process at 
low temperatures is completely independent of the value 
of Jeff and will be present even in the limit Jes/kT —^ oo. 

The main charateristic of a Kasteleyn transition is its 
asymmetry: excitations are only possible at the disor- 
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FIG. 7. Linear magnetic suceptiblity x as a function of tem¬ 
perature for fixed fields (circles). The lines show the same 
magnetic susceptiblity calculated using solely conhgurations 
obeying the spin-ice rule (see text). The inset shows the field 
dependence of the Kasteleyn transition Tk as a function of 
temperature as extracted from Ch and xh- The blue tri¬ 
angles correspond to those determined from the calculation 
of Ch and xh using all states, while the red to those using 
spin-ice configurations only. The solid line is the theoretical 
prediction of Tk{H). 

dered side of the transition. This is shown in Fig. [7] 
where the magnetic susceptiblity is plotted as a func¬ 
tion of temperature at a series of fixed fields as obtained 
from our WLA simulation for a system size of L = 3 (cir¬ 
cles). There, it is clearly seen that at low temperatures, 
while the susceptiblity tends to diverge when Tk is ap¬ 
proached from the disordered side -resembling a second- 
order phase transiton- it is flat on the other side -a be¬ 
haviour more akin to a first-order phase transition. This 
transition was initially termed as “3/2-order” and 
later as “K-type”— . 

A similar asymmetric peak should be expected in the 
specific heat, Ch- Fig. |5]shows Ch as a function of tem¬ 
perature for a series of fixed fields as extracted from our 
simulations. At low fields, it is easy to distinguish two 
clearly defined peaks, the first one, at high temperatures 
corresponds to the Schottky-like peak already mentioned 
in the H — 0 section corresponding to the onset of spin- 
ice correlations in the system. The low temperature peak 
corresponds to the Kasteleyn transition, and shows the 
expected sharp edge at low temperatures and gradual rise 
on the high temperature side. As the field is increased, 
the K-type transition is moved towards higher tempera¬ 
tures. In our simulations Jes/k = l.llK so the condition 
of Jeff "C kT is no longer satisfied at Tk{H) and point like 
excitations are seen at both sides of the transition, grad¬ 
ually changing the peak into a more symmetrical shape. 
As point excitations become more important, the simple 
argument sketched above for the held dependence of Tk 
is no longer valid, and the transition widens and deviates 
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FIG. 8. Specific Heat vs temperature for L = 4 and differ¬ 
ent values of field along [100] (circles). For low fields, it is 
easy to distinguish the high temperature Schottky-type peak 
characteristic of the onset of spin-ice correlations from the 
low temperature asymmetric peak due to the K-type transi¬ 
tion. As the field is raised, the transition moves to higher 
T and it is gradually affected by the presence of additional 
local excitations. The lines show the specific heat calculated 
using only configurations obeying the ice-rule, consequently, 
the Schottky-type peak is absent. 


from a linear dependence (see the blue triangles in the 
inset of Fig. [T]). 

In the WLA, it is possible to impose additional con¬ 
straints on the Hamiltonian by selectively choosing the 
states used to construct the partition function. In this 
case in particular, it is very simple for finite size lat¬ 
tices to identify the states that strictly obey the ice-rule 
by their magnetisation value. In this way, it is possi¬ 
ble to calculate the different thermodynamic quantities 
for the ideal Kasteleyn transition, isolating the effect of 
point-like defects. The results for the susceptiblity and 
specific heat are shown as solid lines in figures [7] and [5] 
respectively. They coincide at low temperatures, when 
kT <C Jeff and the unconstrained curves show a transi¬ 
tion at a lower temperature. The most noticeable change 
is seen in the specific heat, where the Schottky-like peak 
is absent from the ice-rule obeying curves. Furthermore, 
if we repeat the analysis to extract Tk{H) from this lat¬ 
ter set of data, the curve (red triangles in the inset to Fig 
.[71) follows very closely the theoretical prediction (solid 
line). 

One of the unique characteristics of the WLA is that 
it makes it simple to calculate the dependence of the 
free energy as a function of a chosen order parameter 
-in our case the magnetisation. This provides valuable 
information regarding the nature of a phase transition, 
and is particulary interesting for an unusual case such as 
the K-type transition. 

As we have mentioned before, the Kasteleyn transition 
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FIG. 9. The Gibbs potential for the system as a function 
of the magnetisation calculated using the WLA for ice-rule 
configurations at a fixed temperature T = 0.2K and for three 
different fields along [100]: Hi > Hk (blue symbols), H 2 = 
Hk (green symbols), and Hz < Hk (red symbols). The line 
is a guide to the eye. The inset shows the entropy per spin, 
s as a function of the energy for a fixed field of 0.05 T (blue 
dots). 


takes place when J^s/kT is small enough that excitations 
that break the ice-rule are extremely improbable. In this 
case, the energy of the system is a constant, the free 
energy to a purely entropic term and the Gibbs potential 
is given by G = —TS — MH. This resembles the case of a 
simple paramagnet, however, as discussed by Jaubert and 
collaborators (see and a crucial difference arises 
from the ice-rule constraint: if, contrary to the case of the 
paramagnet, this constraint brings the entropy to zero 
at a finite H/kT, this is sufficient to drive a Kasteleyn 
transition in the system. This ad-hoc supposition can be 
put to test using the WLA. The inset of Fig. shows 
the behaviour of the entropy per spin, s as a function 
of the energy in the neighbourhood of s = 0 for a hxed 
field of 0.05 T. As seen in the figure, the slope at which 
the entropy vanishes is indeed finite and, furthermore, 
it is given precisely by l/IV, with Tk the Kasteleyn 
temperature for this field determined from y and C. 

The main panel of the figure shows the Gibbs potential 
as a function of the magnetisation, G{M), at T = 0.2K 
for three different fields along [100]: Hi > Hk, H 2 = Hk 
and Hz < Hk as determined using the WLA for a sys¬ 
tem of L = 4 using only ice-rule configurations. This 
figure captures the characteristic features expected for 
a K-type transition. The low held curve resembles that 
of a paramagnet, with a wide minimum at a non-zero 
magnetisation. As the held is raised towards Hk this 
minimum becomes wider and moves to higher values of 
M,while huctuations increase. At the critical held Hk, 
the system becomes singular, the minimum sits at Mgat 
and the curve becomes hat {dG/dM = 0) in its neigh¬ 
bourhood. For H > Hk, the absolute minimum sits 


at Afsat, the neighbourhood to the minimum is linear, 
with dG/dM hnite and negative, showing the complete 
absence of huctuations in the ordered state. 


IV. CONCLUSIONS 

In this work we have explored by means of the WLA 
the nearest-neighbour spin-ice model, an example of a 
simple classical frustrated model. We have determined 
the value of the residual entropy Sq by doing a hnite size 
analysis of So{L) which can be calculated directly from 
the DoS determined by the WLA for samples of size L^. 
By including the magnetisation as one of the parameters 
of the DoS, we were able to calculate the thermodynamic 
properties of the nnSI model as a function of held. In 
particular, we have investigated the cases where the held 
was applied in the [111] and the [100] directions. In the 
hrst case, we demonstrated that the magnetisation as a 
function of held calculated with the WLA had all the fea¬ 
tures characteristic to this direction, namely, a 1/T rise 
to the kagome-ice plateau at 3.33 /is/Dy ion and a sud¬ 
den jump at H pe Itesla to the fully saturated state. The 
method allowed for a direct calculation of the held depen¬ 
dent entropy, S{H), without the need of any additional 
hxed parameter. As expected, S{H) has a plateau at the 
KI phase with a value that tends to that determined by 
the mapping of the system into dimer conhgurations on 
a honeycomb, and shows a marked peak at the polari¬ 
sation transition. In a similar fashion, we showed that 
when a held is applied along [100], a Kasteleyn transi¬ 
tion takes place between a S' = 0 state and one where 
line-like excitations pierce the system. Additionally, the 
WLA provides information not obtained through other 
methods, such as the free energy as a function of the 
order parameter, which, to our knowledge, has been cal¬ 
culated for the hrst time for a Kasteleyn transition. By 
this mean, we were also able to prove that the ad-hoc 
assumption that the local constraint of the spin-ice rule, 
brings the entropy to zero at a finite H/kT, is perfectly 
valid. We have, furthermore, showed that the WLA al¬ 
lows the computation of the thermodynamic properties 
of the system when additional constraints -not present 
in the Hamiltonian- are put into the system. In particu¬ 
lar, we showed that by selecting only the ice-rule states, 
we could calculate the behaviour of the ideal Kasteleyn 
transition, that is, the one that takes place in the ab¬ 
sence of point defects. In summary, we show that the 
WLA is a very useful tool for simulations of frustrated 
magnetic systems and provides accurate information for 
the thermodynamic properties in equilibrium. As a re¬ 
sult, quantities such as the entropy and the free energy, 
that is cumbersome to obtain through other methods, 
can be easily computed. In addition, the algorithm can 
be used to study the system under the influence of addi¬ 
tional constraints. 
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